/*  sbstempl.h by K.Tsuru  */
// ver. 2.17
/***************************************************************
SN library
Template classes of binary splitting method using recursive
procedure for a series in a form

                      B(0)...B(k-1)
S_L = sum_{k=0}^{L-1} ------------- A_k
                       C(0)...(k)       .

[reference]
1)E.A.Karatsuba
"Fast evaluation of transcendental functions"
Problems of Information Transmission 27, pp.339-360, 1992
Problemy Peredachi Informatsii 27, pp.76-99, 1991 (in Russian)].

2)B.Haible and T.Papaniklaou
"Fast multiprecision evaluation of series of rational numbers"

3)T.Migita, A.Amano, N.Asada and S.Fujino
"Recursive reduction of series for multiple-precision evaluation
and its application to Pi calculation" (in Japanese)
IPSJ Journal Vol.40, pp.4193-4200(1999).
***************************************************************/
#ifndef BINARY_SPLITTING_TEMPLATE_H
#define BINARY_SPLITTING_TEMPLATE_H

template <class BSType> class BinarySplitting {
protected:
    bool isOk;  // can use the values A,B,C // ver 2.18
    BSType A, B, C;
    long L, prec;  // The value of 'L' can be obtained by "ulong L = SCalcInfo::upToTerm;". See below.
public:
    // A virtual function which sets coefficients Ak, Bk and Ck.
    virtual void setABC(long k, BSType& a, BSType& b, BSType& c) = 0;
    // If precision = 0, it does not call roundOff().
    BinarySplitting(long upto, long precision) : isOk(false), L(upto), prec(precision) {
        SCalcInfo::upToTerm = (ulong)L;
    }
    void putTogether();
    virtual ~BinarySplitting(){}
private:
    void recursive(long n0, long n1, BSType &a, BSType &b, BSType &c);
    void roundOff(BSType &a, BSType &b, BSType &c, long f);
public:
    // If "inv == true" return C / A, else A / C.
    SDouble getValue(bool inv = false);
    bool isOK() const { return isOk; }  // ver 2.18
    void getAC(BSType &a, BSType &c) {
        a = A;  c = C;
    }
    BSType getA() const { return A; }
    BSType getB() const { return B; }
    BSType getC() const { return C; }
};

/****************
* Implementation
*****************/

template <class BSType> void BinarySplitting<BSType>::putTogether() {
    recursive(0, L, A, B, C);   isOk = true;
}
template <class BSType> void BinarySplitting<BSType>::
    recursive(long n0, long n1, BSType &a, BSType &b, BSType &c) {

    BSType rA, rB, rC;

    if(n1 - n0 == 1) {  setABC(n0, a, b, c);    return; }

    long midN = (n0 + n1) / 2;

    recursive(n0, midN, a, b, c);       // recursive call left side
    recursive(midN, n1, rA, rB ,rC);    // right side

    a = a * rC + b * rA;    // A(2k) * C(2k+1) + B(2k) * A(2k+1)
    b *= rB;                // B(2k) * B(2k+1)
    c *= rC;                // C(2k) * C(2k+1)

    if(prec) {
        roundOff(a, b, c, prec);    roundOff(rA, rB, rC, prec);
    }
}
template <class BSType> void BinarySplitting<BSType>::
    roundOff(BSType &a, BSType &b, BSType &c, long f) {
/* ver 2.18
    if (a.RawSign() == a.UNDECIDED || b.RawSign() == b.UNDECIDED
        || c.RawSign() == c.UNDECIDED) return;
*/
    long m = (long)a.Head() - f;
    if (m > 0) a.FigureClear(0, (uint)m);
    m = (long)b.Head() - f;
    if (m > 0) b.FigureClear(0, (uint)m);
    m = (long)c.Head() - f;
    if (m > 0) c.FigureClear(0, (uint)m);
}
template <class BSType> SDouble BinarySplitting<BSType>::getValue(bool inv) {
    if (!isOk) putTogether();   // ver 2.18
    if (inv) return SDouble(C) / SDouble(A);
    return SDouble(A) / SDouble(C);
}

/***************************************************************
in the case of B(k) = 1, e.g. e - 1 = 1/1! + 1/2! + 1/3! + ....
See sdcbse.cpp
This can be also used the radix conversion of desimal. Let R radix

S_L =(A_0 . A_1 A_2 A_3...)_R^(-1)
= A_0 + A_1/R + A_2/R^2 + A_3/R^3 + ...
                            A_k
 = sum_{k=0}^{L-1} -------- 
                            R^k    .  
****************************************************************/

template <class BSType> class BinarySplittingA1C {
protected:
    bool isOk;  // can use the values A, C // ver 2.18
    BSType A, C;
    long L, prec;  // The value of 'L' can be obtained by "ulong L = SCalcInfo::upToTerm;". See below.
public:
    // A virtual function which sets coefficients Ak and Ck.
    virtual void setAC(long k, BSType& a, BSType& c) = 0;
    // If precision = 0, it does not call roundOff().
    BinarySplittingA1C(long upto, long precision) : isOk(false), L(upto), prec(precision) {
        SCalcInfo::upToTerm = (ulong)L;
    }
    void putTogether();
    virtual ~BinarySplittingA1C(){}
private:
    void recursive(long n0, long n1, BSType &a, BSType &c);
    void roundOff(BSType &a, BSType &c, long f);
public:
    // If "inv == true" return C / A, else A / C.
    SDouble getValue(bool inv = false);
    bool isOK() const { return isOk; }  // ver 2.18
    void getAC(BSType &a, BSType &c) {
        a = A;  c = C;
    }
    BSType getA() const { return A; }
    BSType getC() const { return C; }
};

/****************
* Implementation
*****************/

template <class BSType> void BinarySplittingA1C<BSType>::putTogether() {
    recursive(0, L, A, C);  isOk = true;
}
template <class BSType> void BinarySplittingA1C<BSType>::
    recursive(long n0, long n1, BSType &a, BSType &c) {

    BSType rA, rC;

    if(n1 - n0 == 1) {  setAC(n0, a, c);    return; }

    long midN = (n0 + n1) / 2;

    recursive(n0, midN, a, c);      // recursive call left side
    recursive(midN, n1, rA, rC);    // right side

    a = a * rC + rA;    // A(2k) * C(2k+1) + A(2k+1)
    c *= rC;            // C(2k) * C(2k+1)

    if(prec) {
        roundOff(a, c, prec);   roundOff(rA, rC, prec);
    }
}
template <class BSType> void BinarySplittingA1C<BSType>::
    roundOff(BSType &a, BSType &c, long f) {

//  if (a.RawSign() == a.UNDECIDED || c.RawSign() == c.UNDECIDED) return; // ver 2.18

    long m = (long)a.Head() - f;
    if (m > 0) a.FigureClear(0, (uint)m);
    m = (long)c.Head() - f;
    if (m > 0) c.FigureClear(0, (uint)m);
}
template <class BSType> SDouble BinarySplittingA1C<BSType>::getValue(bool inv) {
    if (!isOk) putTogether();   // ver 2.18
    if (inv) return SDouble(C) / SDouble(A);
    return SDouble(A) / SDouble(C);
}
/**************************************************************************
in the case of C(k) = 1, e.g. S_L = A_0 + A_1 * B_0 + A_2 * B_0 * B_1 + ...
It can be used to evaluate a polynomial of SLong or SInteger.
**************************************************************************/
template <class BSType> class BinarySplittingPoly {
protected:
    bool isOk;  // can use the values A,B // ver 2.18
    BSType A, B;
    long L; // The value of 'L' can be obtained by "ulong L = SCalcInfo::upToTerm;". See below.
public:
    // A virtual function sets coefficients Ak and Bk.
    virtual void setAB(long k, BSType& a, BSType& b) = 0;
    BinarySplittingPoly(long upto) : isOk(false), L(upto) {
        SCalcInfo::upToTerm = (ulong)L;
    }
    void putTogether();
    virtual ~BinarySplittingPoly(){}
private:
    void recursive(long n0, long n1, BSType &a, BSType &b);
public:
    bool isOK() const { return isOk; }  // ver 2.18
    BSType getValue();  //It returns the value of A.
    BSType getB() const { return B; }
};

/****************
* Implementation
*****************/

template <class BSType> void BinarySplittingPoly<BSType>::putTogether() {
    recursive(0, L, A, B);  isOk = true;
}
template <class BSType> void BinarySplittingPoly<BSType>::
    recursive(long n0, long n1, BSType &a, BSType &b) {

    BSType rA, rB;

    if(n1 - n0 == 1) {  setAB(n0, a, b);    return; }

    long midN = (n0 + n1) / 2;

    recursive(n0, midN, a, b);      // recursive call left side
    recursive(midN, n1, rA, rB);    // right side

    a = a + b * rA; // A(2k) + B(2k) * A(2k+1)
    b *= rB;        // B(2k) * B(2k+1)
}
template <class BSType> BSType BinarySplittingPoly<BSType>::getValue() {
    if (!isOk) putTogether();   // ver 2.18
    return A;
}

/***********************************************************************
It evaluates exp(n/d) with BSType n and d using binary splitting method.
There is a use example in "sdfexbsr.cpp".
************************************************************************/
template <class BSType> class ExpBSRationalNumber
    : public BinarySplitting<BSType> {
    BSType num, den;    // numerator and denominator
public:
    ExpBSRationalNumber(long L, long prec,
        const BSType& n, const BSType& d)
         : BinarySplitting<BSType>(L, prec), num(n), den(d) {}

    void setABC(long k, BSType& a, BSType& b, BSType& c) {
        a.SetInt(1);    b = num;
        if (k == 0) c.SetInt(1);
        else c = k * den;
    }
};

/*********************************************************************
SFraction class is not used. When SFraction class is used, the speed becomes
 slower and the executable file size larger.
**********************************************************************/
struct SLFraction {
    SLong num, den;
};
struct SDFraction {
    SDouble num, den;
};
/*********************************************************************
It splits a SDouble value X in the form of rational number series

 X = u_0 + u_1/v + u_2/v^2 + u_3/v^4 + u_4/v^8 + ... + u_p/v^(2^{p-1}).

where v = S_BASE. It returns 'p'.
Then we obtain
 exp(X) = exp(u_0)exp(u_1/v)exp(u_2/v^2)...
Each term can be evaluated using binariy splitting method.
**********************************************************************/
int SplitSDouble(const SDouble& X, SNBlock <SLFraction>& slr, const SLong& V);  // 3800
#endif // BINARY_SPLITTING_TEMPLATE_H